// M Orsi - May 2007. This program computes the permeation resistance profile and the corresponding permeability coefficient
// Input: a file "GD.dat" of this type:
// [ free energy difference ]  [ diffusion coefficient ]

// Input: file "dz.in", containing the z-distance (in nm) between the selected z-positions 
// Output: resistance function R(z) and permeability P

// Compilation: gcc -Wall -O -o rp rp.c -lm

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#define	kB_J_K		1.3806505e-23 /* [J / K] */
#define AVOGAD		6.0221367e23 /* [mol^(-1)] */
#define kB_J_K_mol      (kB_J_K * AVOGAD) // [ J / (K mol) ]
#define kB_kJ_K_mol      (kB_J_K_mol / 1000) // [ J / (K mol) ]
#define T                303  // temperature [K]
#define kT           ( kB_kJ_K_mol * T ) /* [ kJ / mol) ] */
#define kT2           ( kT * kT ) 

#define AllocMem(a, n, t)  a = (t *) malloc ((n) * sizeof (t))
#define Min( x1, x2 )  ( ( ( x1 ) < ( x2 ) ) ? ( x1 ) : ( x2 ) )

// following defines and application from C: a Reference Manual, Harbison & Steele, 5th ed., p.355; but also Rapaport [2004] p.496.
#define LINELENGTH 80 
#define SEPCHARS " \t\n" 


// global variables
double *deltaG, // free energy array
  *Dz,   // diffusion coefficient
  *R,   // resistance profile
  P;     // permeability coeff

int count, nValues; // number of data-points
double dz;

void AllocArrays()
{
  FILE *fPtr;
  char buff[80];
  // count data
  if ( ( fPtr = fopen( "GD.dat", "r" ) ) == NULL)
    printf( "File rp.dat could not be opened\n" );
  else {
    count = -1;
    while ( !feof( fPtr ) ) {
      fgets( buff, 80, fPtr );
      //printf("tRun = %d\n", tRun);
      count++;
    }  
    fclose( fPtr );
  }
  printf( "count = %d = # data per z = # lines in the rp.dat file\n", count );
  nValues = 2 * count;
  printf( "Number of values for the resistance profile: 2*count = %d, as we symmetrize the data.\n", nValues );
  AllocMem( deltaG, nValues, double );
  AllocMem( Dz, nValues, double );
  AllocMem( R, nValues, double );
  printf( "Allocating memory... done\n" );
}

void  LoadData()
{
  FILE *fPtr;
  char line[ LINELENGTH ], // buffer for each line of the input data file
    *token;        // to store tokens from each "buff line" 
  double d;
  int i;
  
  if ( ( fPtr = fopen( "GD.dat", "r" ) ) == NULL) printf( "File GD.dat could not be opened\n" );
  else {
    i = 0;
    while ( 1 ) {
      fgets( line, LINELENGTH, fPtr );
      if ( feof( fPtr ) ) break;
      if ( strlen( line ) <= 1 ) break;
      token = strtok( line, SEPCHARS ); // find first token (numer in this case)
      d = atof( token );
      deltaG[ i ] = deltaG[ nValues - 1 - i ] = d;
      token = strtok( NULL, SEPCHARS ); // find first token (numer in this case)
      d = atof( token );
      Dz[ i ] = Dz[ nValues - 1 - i ] = d;     
      i++;
    }

    printf( "deltaG[ kJ / mol ]\tDz[ cm^2 / s ]\n" );
    for ( i = 0; i < nValues; i++ ) {
      printf( "%5.3f\t\t%5.3e\n", deltaG[ i ], Dz[ i ] );
    }
  }
}

void  ComputeResistanceProfile()
{
  int i;
  FILE *fPtr;
  
  for ( i = 0; i < nValues; i++ ) {
    R[ i ] = exp( deltaG[ i ] / kT ) / Dz[ i ]; // [ s / cm^2 ]
  }

  // printing 
  if ( ( fPtr = fopen( "rp.dat", "w" ) ) == NULL )
    printf( "File rp.dat could not be opened\n" );
  else { 
    for ( i = 0; i < nValues; i++ ) fprintf( fPtr, "%f\n", R[ i ] ); // [ s / cm^2 ]
    fclose( fPtr );
  }
}

void LoadZstep() 
{
    FILE *fPtr;

    if ( ( fPtr = fopen( "zStep.in", "r" ) ) == NULL) { printf( "File zStep.in could not be opened\n" ); exit( 0 ); }
    else {
      fscanf( fPtr, "%lf", &dz );
      fclose( fPtr );
    }
}
  
void PrintProfiles()
{
  FILE *g, *d, *r;
  int i;

  if ( ( g = fopen( "G.dat", "w" ) ) == NULL ) printf( "File G.dat could not be opened\n" );
  if ( ( d = fopen( "D.dat", "w" ) ) == NULL ) printf( "File D.dat could not be opened\n" );
  if ( ( r = fopen( "R.dat", "w" ) ) == NULL ) printf( "File R.dat could not be opened\n" );

  for ( i = 0; i < nValues; i++ ) {
    fprintf( g, "%f\t%f\n", i * dz - ( dz * nValues / 2 - dz / 2 ), deltaG[ i ] );
    fprintf( d, "%f\t%f\n", i * dz - ( dz * nValues / 2 - dz / 2 ), Dz[ i ] *1E5 );
    fprintf( r, "% f\t%30.20f\n", i * dz - ( dz * nValues / 2 - dz / 2 ), R[ i ] *1E-6 );
  }
  fclose( g );
  fclose( d );
  fclose( r );
}


void IntegrateResistanceProfile() // to get the permeability coefficient
{
  int i;
  double intR, P;
  intR = .5 * ( R[ 0 ] + R[ nValues - 1 ] );
  for ( i = 1; i < nValues - 1; i++ ) intR += R[ i ];
  dz *= 1e-7; // converting dz from nm to cm
  intR *= dz; // space-scaling
  P = 1. / intR;
  printf( "Permeability coefficient P = %.3e cm / s\n", P ); 
  printf( "Permeability coefficient P = %.6f micron / s\n", P * 1e4 ); 
}

int main()
{
  AllocArrays();
  LoadData();
  ComputeResistanceProfile();
  LoadZstep();
  PrintProfiles();
  IntegrateResistanceProfile();
  return 0;
}

